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BACKGROUND OF THE INVENTION 
[0004] The present invention relates to diagnostic ultrasound systems which measure and 
image anatomical structures, and their movement. More particularly, the present invention 
relates to a signal processing method and apparatus for calculation and display of tissue 
deformation to be used in ultrasonic imaging systems. 

[0005] Recently, within the field of ultrasound imaging, physicians have become interested 
in using tissue deformation properties, such as tissue strain and strain velocity for clinical 
measurements. 

[0006] The term "strain" refers to a characteristic of material being examined. For example, 
the strain associated with muscle tissue corresponds to a ratio of the change in muscle tissue 
length during a prescribed time interval to the muscle tissue's initial length. In ultrasound 
imaging, the rate of change of strain (e.g., strain rate, strain velocity, etc.) may be visually 
presented to a physician as a colorized 2-dimensional image, where variations in color 
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correspond to different strain velocities. It has become apparent that the viability of a 
segment of the muscle is related to the amount of muscle strain and temporal behavior of the 
strain that is performed by, or is imposed on the muscle segment. Also, it has been 
determined that malignant tumors may be detected based on their resistance to compression. 

[0007] One application of real-time strain velocity imaging is in cardiology. The strain 
velocity gives a direct and quantitative measure for the ability of the myocardium to contract 
and relax. By imaging along the myocard from an apical view, the local strain velocity 
component along the long axis of the heart can be measured. Measuring the local strain 
velocity component gives information about the local shortening and lengthening of the heart 
wall. By imaging from the parasternal view, the strain velocity component perpendicular to 
the heart wall can be found. Finding the strain velocity component perpendicular to the heart 
wall gives information about the local thickening of the muscle. Wall thickening measured 
with M-mode or from the 2D image is a commonly used measure for muscle viability. With 
strain velocity imaging, a direct measure for this thickening is available. The strain velocity 
images can potentially add to the diagnosis of a number of cardiac disorders. 

[0008] Another application of strain velocity imaging is in heart transplants. Velocity 
variations inside the myocardium are important for the diagnosis of rejection after heart 
transplantation. The strain velocity images give a direct display of these velocity variations. 

[0009] Another application of strain velocity imaging is in noninvasive electrophysiology. 
The preferred embodiment describes techniques to image the local contraction/relaxation 
contributions with a high spatial and temporal resolution. Local contraction/relaxation 
information can be used to accurately determine the localization of, for example, where the 
mechanical movement in the heart chambers is activated based on a cross section just below 
the AV-plane. Furthermore, aberrant conduction pathways (Wolf-Parkinson- White) from the 
atrium to the ventricle can be localized for later ablation. Even the depth inside myocard of 
these paths can be better localized with this invention in order to determine if the patient 
should be treated with catheter techniques or surgical techniques. 

[0010] Another application of strain velocity imaging is in measuring cardiac wall 
thickening. A well established methodology in cardiac diagnosis is to acquire a M-Mode 
image and to measure the wall thickening of myocardium during systole. The preferred 



2 



embodiment provides techniques to take this wall thickening information and measure it in 
real-time with a high precision in both the spatial and temporal domain. The high diagnostic 
relevance of the current wall thickening measurements indicates that the imaging modality 
described in this invention contains highly relevant information for cardiac diagnosis 

[0011] To understand strain velocity or strain rate in more detail, it is assumed that an object 
of initial length Lo may be stretched or compressed or itself lengthens or contracts to a 
different length L. The one-dimensional strain, defined as 

e = L^Jul (i) 

represents a dimensionless description of the change. If the length L is considered to be a 
function of time, the temporal derivative of the strain, the strain velocity, can be found using 
the equation: 

[0012] If the velocity, v of every point in the object is known, an equivalent definition of the 
strain velocity is: 

* = — (3) 
or 

[0013] These equations also provide a useful description of the deformation of the object. In 
Eq. 3, r is the spatial direction of the stretching or compression. The relation between Eq. 2 
and Eq. 3 can be seen if the length L is defined as L(t) = ri(t)-r\(i), and Lo=L(/o), where n is 
the distance to one end of the object, and ri is the distance to the other, and letting t -> tO, and 
letting rl -> r2. As illustrated in Eq. 3, the strain velocity is in fact the spatial gradient of the 
velocity. The strain velocity thus measures the rate of the deformation of the object. If the 
strain velocity is zero, the shape of the object is not changing. If the strain velocity is 
positive, the length of the object is increasing, and if the strain velocity is negative, the length 
of the object is decreasing. Strain velocity is also known as rate-of-deformation, stretching, 
strain rate or velocity strain, 

[0014] Strain imaging is currently an established research area in ultrasonic imaging. The 
degree of deformation in imaged structure may be estimated by correlation of 2D images 
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obtained before and after a pressure increase. One disadvantage of estimating image 
deformation based on correlation of images is that the instantaneous value of the strain is not 
calculated nor displayed in real-time. The lack of a real-time capability is an important 
clinical disadvantage. For example, if strain imaging could be performed in real-time, strain 
imaging could be applied more effectively in cardiac ultrasound or could be used as an 
interactive inspection modality where anomalies in tissue compressibility can be visualized 
in real-time according to the pressure gradients that are applied to the imaged structures. 

[0015] A method of position tracking has been proposed to estimate the local myocardial 
strain velocity based on radio frequency (RF) M-Mode acquisitions. The position tracking 
method is described in H. Kanai, H. Hasegawa, N. Chubachi, Y. Koiwa, and M. Tanaka, 
"Noninvasive evaluation of local myocardial thickening and its color-coded imaging," IEEE 
Trans, on Ultrasonics, Ferroelectrics and Frequency Control, vol. 44, pp. 752-768, 1997. 
However, the method described in the Kanai et al. article has the disadvantages of poor 
temporal resolution and high computational cost, which render real-time imaging difficult 
and costly. Furthermore, the method described in the Kanai et al. article is a manual M-mode 
technique, not well suited to form the basis for real-time two-dimensional strain images. 
Also, the strain velocity is a derivative of a velocity estimate and is therefore very noise 
sensitive. The fundamental velocity aliasing problem that is inherent in tissue velocity 
imaging makes noise difficult to overcome because aliasing prevents the pulse repetition 
frequency from being set at a low enough rate to allow a large observation time. If the 
observation time could be increased, the noise robustness of the strain velocity images could 
be significantly improved. 

[0016] Certain of the above identified difficulties are addressed and overcome according to 
the teachings of United States Patent Application No. 09/167,896, filed October 7, 1998 and 
entitled "A METHOD AND APPARATUS FOR PROVIDING REAL-TIME 
CALCULATION AND DISPLAY OF STRAIN IN ULTRASOUND IMAGING," which is 
incorporated herein by reference. However, is an object of the present invention to 
supplement and/or improve upon such teachings. Certain additional difficulties and 
shortcomings of the prior art are described below. 



4 



[0017] To achieve high frame rate in color Doppler applications, two previously known 
techniques are commonly used: multi line acquisition (MLA) and interleaving. These 
techniques make it possible to acquire more data than in a basic mode, where the scanner 
after one pulse is received waits the specified pulse repetition time (7) before firing the next 
pulse in the same direction. The time to acquire a frame of Doppler data in the basic mode 
is: 



t D0 =N b NT, 



(4) 



where N is the number of pulses in each direction and Nb is the number of beams in the 
image. A relatively small extra delay related to the change in setup of the transmitter and 
beamformer is ignored to simplify the discussion. 

[0018] In the MLA method, a broad beam is transmitted. When receiving the echo, the 
signals from all the transducer elements are processed in parallel in two or more 
beamformers. Each beamformer time delays the element signals differently to generate 
different receive beams. This way, two or more beams can be acquired during the time for 
one pulse-echo cycle, and the frame rate can be increased correspondingly. Using MLA, the 
time to acquire a frame of Doppler data is 

t DM LA=p-NT, (5) 

" MLA 

where N M u is the number of beams that are processed in parallel. 

[0019] In the interleaving technique, the waiting time T from one pulse to the next in the 
same direction is utilized to send pulses in other directions, as illustrated in Fig. 1. There is 
however a minimum waiting time T 0 where no other pulses can be fired in any direction. 
This is given by the time for the pulse to travel to the maximum depth and back: T 0 > 2d/c. 
The number of directions that pulses are fired during the time T is called the interleave group 
size, N int . This obviously has to be an integer number, and T = N int T 0 . Using interleaving, 
the time to acquire a frame of Doppler data becomes: 
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Dint 




NT. 



(6) 



[0020] Fig. 1 illustrates the pulse order and beam directions in the interleaving method for 
three different interleave group sizes N inh The number of beams, Nb, equals 8 and packet 
size, N, equals 2 in the example of Fig. 1. For interleave pattern 100, the interleave group 
size, N inh equals 8, for interleave pattern 110, N int equals 4 and for interleave pattern 120, N int 
equals 1. 

[0021] A typical scanning procedure for a tissue Doppler application is illustrated in Fig. 2. 
In the example of Fig. 2, the packet size N equals 3 and the interleave group size Nint equals 
Nb. T is the pulse repetition time, tj and to are the times needed to acquire a tissue frame and 
a Doppler frame respectively, and t F is the total acquisition time for one tissue Doppler 
frame. A tissue frame 130 is first captured, using high beam density. The PRF used for 
tissue Doppler is usually so low that only one interleave group is necessary. So N Doppler 
sub frames 132, 134 and 136 are captured separately, usually using fewer beams than in the 
tissue frame. The velocity is calculated from the N sub frames 132, 134 and 136, color coded 
and then mapped onto the tissue frame. The time to acquire a tissue Doppler frame then 
becomes: 



where tj is the time required to acquire the tissue frame. It is, thus, apparent that the 
maximum frame rate is limited by the above described ultrasound data acquisition schemes. 

[0022] It is known that tissue velocity can be estimated using either the first or the second 
harmonic component of the ultrasound signal. Using the second harmonic component 
{octave imaging) has been reported to produce an improvement in image quality in gray scale 
images, and the same improvement can be expected in tissue Doppler. There is however a 
disadvantage in that the Nyquist limit is halved when using the second harmonic instead of 
the fundamental component. Using a low PRF is also preferable, since the phase amplitude 
of the complex signal is increased compared to the noise, resulting in a lower variance in the 
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velocity estimate. A disadvantage of using a low PRF is that the Nyquist limit is further 
reduced. A reduced Nyquist limit increases the risk of aliasing, which results in the 
misrepresentation of high velocities. 

BRIEF SUMMARY OF THE INVENTION 
[0023] An ultrasound system and method for calculation and display of tissue deformation 
parameters are disclosed. 

[0024] According to one aspect of a preferred embodiment of the present invention, an 
ultrasound acquisition technique that allows a high frame rate in tissue velocity imaging or 
strain rate imaging is disclosed. With this acquisition technique the same ultrasound pulses 
are used for the tissue image and the Doppler based image. A sliding window technique is 
used for processing. 

[0025] According to another aspect of a preferred embodiment of the present invention, 
strain is determined by an accumulation of strain rate estimates for consecutive frames over 
an interval. The interval may be a triggered interval generated by, for example, an R-wave in 
an ECG trace. The strain calculation may be improved by moving the sample volume from 
which the strain rate is accumulated from frame-to-frame according to relative displacement 
of the tissue within the original sample volume. The relative displacement of the tissue is 
determined by the instantaneous tissue velocity of the sample volume. 

[0026] According to another aspect of a preferred embodiment of the present invention, dr, 
the spatial offset used in an estimation of strain rate, is varied adaptively throughout the 
image. The spatial offset, dr, can be maximized to cover the entire tissue segment (e.g., heart 
wall width) while still keeping both of the sample volumes at each end of the offset within 
the tissue segment. This may be accomplished by determining whether various parameters 
(e.g., grayscale value, absolute power estimate, magnitude of the autocorrelation function 
with unity temporal lag and/or magnitude of strain correlation) of the sample volumes within 
in the spatial offset are above a given threshold. 

[0027] According to another aspect of a preferred embodiment of the present invention, a 
generalized strain rate estimator that is based on a weighted sum of two-sample strain rate 
estimators with different spatial offsets is employed. The weights are proportional to the 
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magnitude of the strain rate correlation estimate for each spatial offset, and thus reduce the 
effect of noisy, i.e. poorly correlated, samples. 

[0028] According to another aspect of a preferred embodiment of the present invention, an 
improved signal correlation estimator that uses a spatial lag in addition to the usual temporal 
lag is disclosed. The spatial lag is found from the tissue velocity. The improved signal 
correlation estimator can be utilized both in the estimation of strain rate and tissue velocity. 

[0029] According to another aspect of a preferred embodiment of the present invention, 
tissue velocity is estimated in a manner that reduces aliasing while maintaining spatial 
resolution. Three copies of a received ultrasound signal are bandpass filtered at three center 
frequencies. The middle of the three center frequencies is centered at the second harmonic of 
the ultrasound signal. A reference tissue velocity is estimated from the two signals filtered at 
the outside center frequencies. The reference tissue velocity is used to choose a tissue 
velocity from a number of tissue velocities estimated from the signal centered at the second 
harmonic. 

[0030] According to another aspect of a preferred embodiment of the present invention, a 
method to estimate the strain rate in any direction, not necessarily along the ultrasound beam, 
based on tissue velocity data from a small region of interest around a sample volume is 
disclosed. 

[0031] According to another aspect of a preferred embodiment of the present invention, a 
plurality of quantitative tissue deformation parameters, such as tissue velocity, tissue velocity 
integrals, strain rate and/or strain, may be presented as functions of time and/or spatial 
position for applications such as stress echo. For example, strain rate or strain values for 
three different stress levels may be plotted together with respect to time over a cardiac cycle. 
Parameters which are derived from strain rate or strain velocity, such as peak systolic wall 
thickening percentage, may be plotted with respect to various stress levels. 

[0032] Other objects, features, and advantages of the present invention will be apparent from 
the accompanying drawings and from the detailed description that follows below. 
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BRIEF DESCRIPTION OF SEVERAL VIEWS OF THE DRAWINGS 
[0033] Fig. 1 illustrates the pulse order and beam direction for three different interleave 
group sizes. 

[0034] Fig. 2 illustrates a typical ultrasound acquisition procedure for a tissue Doppler 
application. 

[0035] Fig. 3 illustrates an ultrasound system according to a preferred embodiment of the 
present invention. 

[0036] Fig. 4 illustrates an ultrasound acquisition procedure for a tissue Doppler, strain or 
strain rate application according to a preferred embodiment of the present invention. 

[0037] Fig. 5 illustrates a graph of the results of a computer simulated comparison of a linear 
regression strain rate estimator according to the prior art and a weighted strain rate estimator 
according to a preferred embodiment of the present invention. 

[0038] Fig. 6 illustrates a graph of the results of a computer simulated comparison of a linear 
regression strain rate estimator according to the prior art and a weighted strain rate estimator 
according to a preferred embodiment of the present invention. 

[0039] Fig. 7 illustrates the coordinates r, u, v, and w, and the insonation angle a used by a 
strain rate estimate angle correction technique according to a preferred embodiment of the 
present invention. 

[0040] Fig. 8 illustrates the velocity components vv, vw and vr, the distance Dr and the angle 
a in a small muscle segment used by a strain rate estimate angle correction technique 
according to a preferred embodiment of the present invention. 

[0041] Fig. 9 illustrates a display of strain rate from multiple stress levels as a function of 
time for a normal case according to a preferred embodiment of the present invention. 

[0042] Fig. 10 illustrates a display of accumulated strain from multiple stress levels as a 
function of time for a normal case according to a preferred embodiment of the present 
invention. 

[0043] Fig. 1 1 illustrates a display of strain rate from multiple stress levels as a function of 
time for a pathological case according to a preferred embodiment of the present invention. 
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[0044] Fig. 12 illustrates a display of accumulated strain from multiple stress levels as a 
function of time for a pathological case according to a preferred embodiment of the present 
invention. 

[0045] Fig. 13 illustrates a display of a strain derived parameter, peak systolic wall 
thickening, as a function of stress level according to a preferred embodiment of the present 
invention. 

DETAILED DESCRIPTION OF THE INVENTION 
[0046] A method and apparatus are described for generating diagnostic images of tissue 
deformation parameters, such as strain rate, strain and tissue velocity, in real time and/or in a 
post-processing mode. In the following description, numerous specific details are set forth in 
order to provide a thorough understanding of the preferred embodiments of the present 
invention. It will be apparent, however, to one of ordinary skill in the art that the present 
invention may be practiced without these specific details. 

[0047] A block diagram for an ultrasound imaging system according to a preferred 
embodiment of the present invention is shown in Fig. 3. A transmitter 140 drives an 
ultrasonic transducer 142 to emit a pulsed ultrasonic beam 144 into the body. The ultrasonic 
pulses are backscattered from structures in the body, like muscular tissue, to produce echoes 
which return to and are detected by the transducer 142. A receiver 146 detects the echoes. 
The echoes are passed from the receiver 146 to a complex demodulation stage 148 and a 
tissue processing stage 149. The complex demodulation stage 148 demodulates the echo 
signals to form I, Q data pairs representative of echo signals. The demodulated I, Q data 
pairs are complex Doppler signals that are passed to a tissue deformation calculation stage 
150 which carries out tissue velocity, strain rate and/or strain calculations as explained 
below. The complex Doppler signal is associated with a sample volume defined by a range 
position and beam in a region of interest. A complex Doppler signal typically comprises a 
segment of data samples which is used to estimate the Doppler shift. The echo signals are 
also passed to the tissue processing stage 149, which performs processing such as B-mode 
processing to form a 2D or 3D image of the anatomical structure scanned. 

[0048] The tissue deformation values, e.g., tissue velocity, strain rate and/or strain, output by 
the tissue deformation calculation stage 150 and the tissue image values output by the tissue 
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processing stage 149 are passed to a display system 152 for display. The display system 152 
. includes a monitor 1 54. 

[0049] United States Patent Application No. 09/167,896, filed October 7, 1998 and entitled 
"A METHOD AND APPARATUS FOR PROVIDING REAL-TIME CALCULATION 
AND DISPLAY OF STRAIN IN ULTRASOUND IMAGING/ 5 which is incorporated herein 
by reference, describes a manner in which a strain rate may be estimated using the system of 
Fig. 3. 

[0050] For Strain Rate Imaging (SRI) and other Doppler based applications where a low 
pulse repetition frequency (PRF) is acceptable, a scanning procedure that allows higher 
frame rate may be used. Instead of collecting separate tissue frames as illustrated in Fig. 2, 
the number of beams in the Doppler sub frames can be increased to allow tissue visualization 
based on only these frames. The acquisition of separate tissue frames becomes unnecessary. 
Fig. 4 illustrates a scanning procedure that allows a high frame rate. This scanning procedure 
may be used in either tissue Doppler or SRI applications. In the example of Fig. 4, the 
packet size N = 3 and the interleave group size Nj nt = Nb. T is the pulse repetition time, tj 
and to are the times needed to acquire a tissue frame and a Doppler frame respectively, and tp 
is the total acquisition time for one tissue Doppler or SRI frame. As illustrated in Fig. 4, a 
Doppler frame is still generated from N sub frames (the subframes are numbered 160, 161, 
162, 163 and 164), but a sliding window technique may be used, so the time to produce one 
Doppler or SRI frame will be only 

^FSRJ = *T » ( 8 ) 

assuming that the time to acquire one Doppler subframe is equal to the time to acquire one 
tissue frame in the conventional method. Comparing equations (7) and (8) one can see that 
the acquisition time for one frame is greatly reduced and, thus, allowing a higher frame rate. 

[0051] One parameter that may be calculated by the tissue deformation calculation stage 150 
is strain. The relation between the strain and the strain rate can be developed by way of an 
example. Consider a one-dimensional homogeneous object of length L(t) that has a spatially 
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constant strain rate field s(t). The term "strain rate" is here used for the velocity gradient. 
The velocity field is thus given as: 

v(t 9 r) = s(t)r, 

(9) 

where r is the position in the object. The velocity at r = 0 is set to zero for simplicity, but the 
same relations will apply also when v(t,0) differs from zero. 

[0052] The change in length over a small time step At can then be estimated as 

L{t + AO - L(t) * Ats{t)L{t). 

(10) 

Letting At —> 0 we get the temporal derivative of the length: 
dL{t) v L(t + At)-L(t) , w/x 

— = hm f — = s(t)L(t). ( li ) 

dt At v ' 

The solution to this differential equation is 

L(0 = 4exp(f; o 5(r)rfr) (12) 

and the strain is finally found as 

e(t) = m ~ L ° -100% = [exp([ o s(r)dr)- 1\ 100%. ( U ) 

[0053] The strain e(i) in a sample volume in the image can be estimated in real-time by 
replacing the integration in equation (13) with a cumulative sum: 

e(0 = [exp(C(0)-l]l00%, 

C(0 = C(/-l) + j(/)Af. K } 

Here i is the frame number and At is the time between each frame. C(i) is the cumulative 
sum, and s(i) is the strain rate estimate for the given sample volume. The accumulation can 
also be reset at any time, for instance at a specific time trigged by an ECG-signal, by setting 
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C(i-l) to zero for the corresponding frame number /. The calculation above can be performed 
for every sample volume in the image, and the visualization can be performed in the same 
way as for tissue velocity imageing (TVI) and SRI, only using a color map representing 
strain rather than tissue velocity or strain rate. 

[0054] A further improvement is possible if the tissue velocity v is also available for each 
sample volume. In the cumulative sum for radial sample volume number m 0 , the strain rate 
estimate might then be taken from a different sample volume given by the tissue velocity. 
First, the frame-to-frame relative displacement index is estimated as 

d = vAtk s ( 15 ) 

where v is the tissue velocity estimate in sample number mo, and k s is the spatial sampling 
frequency. Next, the strain rate estimate from the sample volume number 

m = m Q +d (16) 

is used in the cumulative sum, rather than mo. If the tissue movement is only in the direction 
of the beam, this method allows the cumulative summation to track the motion of the same 
anatomical sample during its movement. Even if the tissue movement is in other directions, 
an improvement is expected. 

[0055] The strain rate estimator in Application No. 09/167,896 was in its simplest form 
described as: 

s(r) = (v(r + dr)-v(r) ) I dr (i 7 ) 

where r is the radial position along an ultrasound beam, v is the tissue velocity, and dr is the 
spatial offset. This spatial offset can be varied adaptively throughout the image. Given an 
upper and lower limit on the size of dr, it can be increased as much as possible while still 
keeping both of the sample volumes at each end of the offset within the tissue. There are 
several different criteria that can be used to ensure that the offset is within the tissue. One 
possible criteria is that the corresponding tissue sample volumes must have a grayscale value 
above a given limit. Another possible criteria is that the power estimates of the sample 
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volumes must have absolute values above a given limit. Another possible criteria is that in 
either of the two sample volumes, the magnitude of the autocorrelation function with unity 
temporal lag must be above a given limit. Another possible criteria is that the magnitude of 
the strain correlation (described in equation (8) in Application No. 09/167,896) must be 
above a given limit. Any of these criteria one can be used separately, or they can be 
combined so that two or more criteria must be met for a positive determination that the 
sample volumes at the end of the offset dr are within the tissue. 

[0056] The tissue deformation calculation stage 14 may calculate strain rate using a strain 
rate estimator that is based on several samples, and is weighted with the magnitude of a strain 
correlation estimate. Consider a quadrature demodulated Doppler signal x(m,n), where m is 
the spatial sample volume index, and n is the temporal index. The signal is assumed to have 
been acquired using a center frequency a pulse repetition time 7, and a radial sampling 
frequency r s equal to the radial size of the point spread function. The speed of sound in the 
imaged object is assumed to be c. An estimator for strain rate based on M spatial and N 
temporal samples of the Doppler signal is: 
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(19) 



is the strain rate correlation estimate, 



d>»=-Z§(m), 



(20) 
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is the angle of the strain rate correlation estimate, and 
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is a weighing factor. The signal correlation estimate &(m) is described below. 

[0057] The strain rate estimator of equation (18) has certain advantages over the prior art 
Myocardial Velocity Gradient (MVG) technique first described in D. Fleming et al., 
"Myocardial velocity gradients detected by Doppler imaging" Br. J. Radiol., 67(799):679- 
688, 1994, and further developed by Uematsu et al., "Myocardial velocity gradient as a new 
indicator of regional left ventricular contraction: Detection by a two-dimensional tissue 
Doppler imaging technique" J. Am. Col. Cardiol., 26(l):217-23, 1995. For example, 
Fleming and Uematsu disclose the use of a least squares linear regression of the velocity data 
to obtain the velocity gradient (strain rate). Linear regression puts equal weight to all the 
velocity samples. The weighted strain estimator of equation (18), however, uses weights that 
vary with the magnitude of the strain rate correlation of equation (19) resulting in an 
improved strain rate estimation. 

[0058] Figs. 5 and 6 illustrate a computer simulation comparison of a least squares linear 
regression estimator and the strain rate estimator of equation (18). Fig. 5 illustrates a linear 
regression fit (dashed line) and a weighted strain rate linear fit (solid line) for simulated 
velocity estimates (circles) at varying depths. Signals including noise were generated with a 
velocity gradient (strain rate) of 1.0 s' 1 . A typical outcome is presented in Fig. 5. Note that 
the two outermost points give a large error for the linear regression line (dashed line), while 
the effect on the weighted strain rate estimator is much less. In Fig. 6, strain rates estimated 
with the linear regression method (stars) and with the weighted strain weight estimator 
(circles) are compared for 50 independent simulations. Once again, signals including noise 
were generated with a velocity gradient (strain rate) of 1.0 s" 1 . The weighted strain rate 
estimator shows less variance than the linear regression method. 

[0059] The signal correlation k{m) (used in equation (19) above) can be estimated in 
different ways. For example, one estimate is 

N-\ 

M™) = JV(m,«)x(m,/i + 1). ( 22 ) 

Spatial averaging may also be used to reduce the variance of k{m) in equation (22) and 
other estimators of k{m) described herein. 
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[0060] A more robust method to estimate the signal correlation &(m) is to introduce a spatial 
lag Am , and correlate signal samples from not just the same depth m, but also from m +Am : 



AM 



R(m) = £ x* (m, ri)x(m + Am, n + 1). ( 23 ) 



n=] 



[0061] The spatial lag Am preferably is chosen to maximize the magnitude of R(m) . One 
way to chose a Am is through a phase root seeking technique such as described in A. 
Peasvento and H. Ermert, "Time-efficient and exact algorithms for adaptive temporal 
stretching and 2D-correlation for elastographic imaging using phase information" Proc. of 
the 1998 Ultrasonic Symposium, to be published, 1998. Alternatively, the inventors have 
found that the peak magnitude of k{m) is found when the spatial lag Am is chosen equal to 
the translation of the tissue from pulse to pulse: 

where v is the tissue velocity, PRF is the pulse repetition frequency and k s is the spatial 
sampling frequency of the signal. This method requires that an un-aliased velocity estimate 
is available. 

[0062] The tissue deformation calculation stage 150 may calculate a velocity estimate as 
follows. Three equal copies of the received signal are band pass filtered with three different 
filters. Two narrow band filters centered at // and y} 3 and a third wider band filter centered at 
fs are used, where fi<fs < fi, and fs is centered around the second harmonic component of the 
signal. The signal correlation of each of these three signals are estimated using equation ( 22 
), resulting in the correlation estimates , k 2 {m) and ft 3 (m) , respectively. The tissue 

velocity can be found from the angle of $ 3 (m) as: 

. cPRF 6 

"' = i^r^ !(m) - < 2s > 

where c is the speed of sound. Unfortunately, the velocity estimate of equation (25) is easily 
aliased. A difference correlation is next found as 
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A d (m) = A l \m)R 2 (m). 



(26) 



[0063] The velocity of the tissue is found from the angle of this difference correlation as 

cPRF h 

rf= ^ra d{mX (27) 

where c is the speed' of sound. This velocity estimate is not as easily aliased since (/}-//)< 
fs. However, with equation (27) the spatial resolution is poor since narrow band signals were 
used in the estimation of k x {m) and k 2 (m). To this point, this two-frequency velocity 
estimation method is similar to a method described in Dousse et al., "Two years experience 
in measuring velocities beyond the Nyquist limit with Color Flow Mapper" Proceedings of 
EURODOP'92, page 219, Brighton, United Kingdom, 1992. 

[0064] To regain the spatial resolution of the estimate in equation (25), the following 
algorithm is used: For each (possibly aliased) velocity estimate i> 3 , several candidate 
velocities are found as 



^=^(^("0+2*4 - 



2{f 2 -A) 



< k < 



A-(f2-A) 

I 2(f 2 -f t ) J 



(28) 



[0065] Next, the candidate velocity \ k that is closest to the (unaliased) difference velocity 
estimate Q d is chosen as the output velocity estimate. This way, the spatial resolution of the 
i> 3 estimates is kept, while avoiding the problem of aliasing. 

[0066] A method for angle correction of a strain rate estimation is described with respect to 
Fig. 7. Locally for each muscle segment of the left cardiac ventricle, the coordinates are 
defined as: 

r - along the ultrasound beam, positive away from the transducer 

1 - lateral (beam-to-beam), positive from left to right in the ultrasound image 

u - circumferential, clockwise seen from the apex 

v - meridional (longitudinal), from apex to base 
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w - transmural, from endo- to epi-card 

where u, v and w will be approximately perpendicular, as shown in Fig. 7. The strain rates in 
these directions are termed s r , s\, s U9 s v and s w , respectively. The origin (u,v,w)=(0,0,0) does 
not need to be defined in relation to the macroscopic ventricle geometry, and can be chosen 
anywhere in the imaged muscle segment. 

[0067] Furthermore, a is defined as the angle between the v-axis and the r-axis, so that zero 
degrees corresponds to measuring along the muscle in the meridional direction. It is assumed 
that the angle a is in the v-w-plane (long axis or apical views), so the problem becomes two- 
dimensional. Note that the angle a is negative in Fig. 7. 

[0068] Without loosing generality, it can be assumed that the point (v,w)=(0,0) is not 
moving. If the strain rate is spatially homogeneous over a small distance Ar in the muscle 
segment, the muscle point (v,w) will then move with the velocity components: 

v v = V5 V ( 29 ) 

and 

v w=w w . (30) 

[0069] These velocity components are shown in Fig. 8. Fig. 8 is an illustration of the 
velocity components v v , v w and v r , the distance Ar and the angle a in a small muscle segment. 
All the parameters are drawn positive, but notice that the angle a is usually negative when 
imaging from the apex, and that v v , and consequently v r , normally are negative during 
systole. The lateral (beam-to-beam) 1-axis is also included for reference. The velocity 
component along the ultrasound beam in position (v,w) becomes 

v r = V5 V cosa + ws w sina. ( 31 ) 
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[0070] Notice that the velocity v r for simplicity is defined positive away from the transducer, 
i.e., in positive r-direction. This is opposite of the usual definition for the velocity sign in 
Doppler imaging. 

[0071] By using velocity-information from more than one beam at a time, it is possible to 
calculate the strain rate in other directions than along the beam. The beams are assumed to be 
parallel in the region of interest. The vw-axis system is then a rotation of the lr-axis system 
by an angle of (a-7i/2), and one can write 

v = rcos# + lsina 

(32) 

w = rsina -lcosa 
[0072] Substituting these equations in equation (31) one gets 

v r = 5 v (rcosor + lsina)cosQ:H-5 w (rsina-lcosa)sina (33) 

[0073] Taking the derivatives in the two directions r and 1, one gets the two equations 

2 2 

— L = s v cos a + s w sin a 

dr (34) 

dv T 

— - = s v sinacosa-s w sinacosa 
51 



Solving for s v and s w gives 



dv r dv T 
s v = — - + — L tana 



dr dl ( 35 ) 



s w = — L -cota 

w dr d\ 

[0074] This means the strain rates in the anatomical directions v (meridional) and w 
(transmural) can be found from the radial and lateral gradients of the measured radial 
velocity, as long as the angle a is known. The image plane lr must coincide with the vw 
plane, which is the case for all apical views and the parasternal long axis view (PLAX). 
Notice that when imaging from the apex, the angle a will be close to zero for most of the 
ventricle. 

[0075] The same formulas apply if one substitutes v with u, so the strain rate in the u- 
direction (circumferential) can also be found. The image plane lr must then coincide with the 
uw plane, which is the case for the short axis view (SAX). 
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[0076] There will be some angles where the strain rates are unavailable, though. For the u or 
the v directions these are the angles where tan a approaches infinite values. For the w- 
direction these are the angles where cot a approaches infinite values. 

[0077] In the SAX view and using a sector scan, an approximation of a can automatically be 
found if the user defines the center of the ventricle. By assuming that the SAX cross section 
of the ventricle is circular, a at a particular location is given as 

<* = Y~ b * c (36) 

where fy, is the angle of the ultrasound beam that intersects the point {Ob = 0 is defined as the 
center beam), and 6 C is the angle between the center ultrasound beam an imaginary beam 
from the center of the ventricle through the point. 

[0078] A preliminary test has been performed using this two-dimensional angle correction 
method. A velocity data set from a healthy volunteer was obtained using tissue Doppler 
imaging with high beam density. The short axis view was used and the circumferential and 
transmural strain rate components in three phases of the cardiac cycle (mid systole, early 
diastolic relaxation and mid diastole) were estimated. The myocardium was segmented 
manually. As expected, the resulting images showed that the radial strain rate is equal to the 
transmural strain rate at twelve and six o'clock, and the circumferential strain rate at two and 
ten o'clock. Except from where the cot a or tan a approach infinity, the apparent noise in the 
images did not seem to be increased by this procedure. 

[0079] It is also possible to perform three-dimensional angle correction. Locally for each 
muscle segment of the left cardiac ventricle, the coordinates are defined as: 

x - azimuthal (perpendicular to the image plane) 

y - lateral (beam-to-beam) 

z - along the ultrasound beam 

u - circumferential, clockwise seen from the apex 

v - meridional (longitudinal), from apex to base 
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w - transmural, from endo- to epi-card 

where the triplets x, y, z and u, v, w locally are assumed to be perpendicular. The strain rates 
in these directions are termed s u , s v and s w , respectively. The origin (u,v,w)=(0,0,0) does not 
need to be defined in relation to the macroscopic ventricle geometry, and can be chosen 
anywhere in the imaged muscle segment. 

[0080] Without loosing generality, it is assume that the point (u,v,w)=(0,0,0) is not moving. 
If the strain rate is spatially homogeneous over a small distance Ar in the muscle segment, 
the muscle point (u,v,w) will then move with the velocity components: 

V u =U *u> V v=™v> V w= W *W (37) 

[0081] Using velocity-information from more than one beam at the time it is possible to 
calculate the strain rate in other directions than along the beam. The beams are assumed to be 
parallel in the region of interest. 

[0082] Based on formulas for axis rotation it is possible to express the velocity components 
in the xyz-directions rather than in the uvw-directions. The velocity component in the z- 
direction (along the ultrasound beam), v z , is the one found using Tissue Velocity Imaging. 
The gradients of this velocity component in each of the three spatial directions are 

v *r=-^:> r = x,y,z (38) 
The relation to the strain rates in the uvw-directions is 



where A(a,/3,y) is a matrix that describes the 3D axis rotation between the uvw-system and 

the xyz-system, and a , J3 , and y are the rotation angles about the u-, v- and w-axis 

respectively. Except for certain rotation angles, this matrix can be inverted, and the strain 
rates can be found as: 
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(40) 



[0083] Estimates for the strain rates in the uvw-directions are found by inserting velocity 
gradient estimators based on recorded tissue velocity data. Examples of estimators for the 
velocity gradients are 



V zr = 



v 2 (r + Ar)-v 2 (r) 
Ar 



r = x,y,z 



(41) 



where Ax, Ay, and Az are the sampling distances in the azimuth, lateral and radial 
directions in the ultrasound data respectively. Similar methods as described for ID strain 
rate can also be used to estimate these velocity gradients, where the radial increment is 
replaced by an increment in the x- and y-directions. 

[0084] A further improvement can be achieved by performing a least squares fit of the 
estimated strain rates to the incompressibility equation 



s t , + s„ + s„, = 0 



(42) 



since muscular tissue can be considered incompressible. 
[0085] In two dimensions the strain rate estimates reduce to 

*u =v 2Z +v^cot£ 
*w =v 2Z tan^ 

for images in the uw-plane (short axis images) and 



S v =V zz + V zy COta 
*w = V zz + V zy tail « 



(43) 



(44) 



for images in the vw-plane (apical images). There will be some angles where the strain rates 
are unavailable, though. For the u or the v directions these are the angles where tan 
approaches infinite values. For the w-direction these are the angles where cot approaches 
infinite values. 
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[0086] The tissue deformation calculations described herein are suited for quantitative stress 
echo applications. There are at least four main quantitative parameters that may be extracted, 
including: tissue velocity, which quantifies wall motion; tissue velocity time integrals which 
quantify accumulated wall motion during a time interval such as systole; strain rate (velocity 
gradient), which quantifies the local wall thickening at a given time instant; and strain 
(integrated strain rate) which quantifies local wall thickening during a time interval such as 
systole. These parameters are functions of both spatial position and time. From these 
parameters, other clinically relevant parameters may be derived. One way to present these 
parameters is to plot pairs of the parameters against each other (similar to pressure-volume- 
loops). Another useful way to present these parameters is to estimate and record (in a 
cineloop for example) one or more parameters from different stress levels in the stress test 
and then display the respective parameters from the varying stress levels simultaneously. 

[0087] During a stress echo examination one of the crucial things to assess is segmental wall 
motion. Typically, the left ventricle is subdivided into segmental areas and a visual 
assessment of wall motion is done in each of these segments from the various cineloops that 
are acquired. The 16 segment ASE model of the left ventricle is currently the most common 
way to subdivide the left ventricle for scoring of stress echo exams. In a visual assessment a 
given segment is compared in terms of motion and wall thickening by visual comparison of 
similar views (2-chamber, 4-chamber, lax, sax, etc.) at different stress levels. The stress 
levels usually include rest, 1 or more levels of stress induced by exercise or pharmacological 
infusion and finally recovery. A normal reading of a segment is that both wall motion and 
local wall thickening increase during systole as a function of the applied stress level. 

[0088] Fig. 9 illustrates how time traces of strain rate for a given location or wall segment 
can be combined from multiple stress levels. Strain rates estimated during rest (line 200), 
medium stress (line 202) and peak stress (line 204) are plotted with respect to time. An ECG 
trace (line 206) is provided for reference at the bottom of the display. A difference in heart 
rate from the various stress levels is accounted for in this example by time scaling the 
different strain rate traces. This combined display contains useful clinical information about 
the local wall function and how the wall segment responds to an increase in stress level. The 
example is a typical normal reading of longitudinal shortening that can be recorded with an 
apical view. It should be noted that the longitudinal shortening assessed in this manner also 
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indirectly describes wall thickening in short axis views because of conservation of mass and 
incompressibility of myocard. The example illustrates a normal reading where longitudinal 
shortening increases with the stress level. 

[0089] Fig. 10 is identical to Fig. 9 except that accumulated strain is plotted for rest (line 
210), medium stress (line 212) and peak stress (line 214) instead of the instantaneous strain 
rate. The Fig. 10 demonstrates how the longitudinal shortening increases as a function of 
stress level. Figs. 11 and 12 correspond to Figs. 9 and 10, respectively, except that Fig. 11 
illustrates a typical pathological reading of strain rates for rest (line 230), medium stress (line 
232) and peak stress (line 234), and Fig. 12 illustrates a typical pathological reading of 
accumulated strain for rest (line 240), medium stress (line 242) and peak stress (line 244). 
The example of Figs. 11 and 12 illustrates a case with normal rest values for longitudinal 
shortening, but with a reduction in shortening when the stress level increases. At peak stress 
the curves illustrate a reverse in both strain rate and strain which can indicate passive 
stretching of the local wall segment. 

[0090] Fig. 13 illustrates how characteristic values extracted from the strain information for a 
given location or wall segment can be displayed as a function of stress level. The example in 
Fig. 13 is the maximal systolic longitudinal shortening which is plotted as a function of stress 
level. The normal case (line 250) with a uniform increase in longitudinal shortening and the 
pathological case (line 252) with a decrease in longitudinal shortening and even a switch to 
passive stretching during systole are illustrated. 

[0091] Another useful way to present the quantitative parameters derived from strain is in a 
Bulls-eye plot by either numerically or graphically labeling each of the areas corresponding 
to LV segments according to the associated strain derived parameters. The values illustrated 
in Fig. 13 are examples of such useful strain derived parameters. 

[0092] In the foregoing specification the invention has been described with reference to 
specific exemplary embodiments thereof. It will, however, be evident that various 
modifications and changes may be made thereto without departing from the broader spirit 
and scope of the invention as set forth in the appended claims. The specification and 
drawings are, accordingly, to be regarding in an illustrative rather than restrictive sense. 
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